*
clear
use "paper_df2.dta"
tsset prio_grid paneldate
***************************************************************************************************
***** CONTROL VARIABLES *****
*****************************************************************************************************
gen poplog = log(pop + 1)
gen gdppclog = log(gdppc + 1)
global X1 pop_log_ mnt_ ttime_log_  nlights_mean
global X2 polity2 gdppclog poplog

gen units_deployedsq = units_deployed^2
gen countries_deployedsq = countries_deployed^2
gen pko_deployedsq = pko_deployed^2
gen lsum_notroopss1 = lsum_notroops^2

gen death_totallog = log(death_total)
gen duration_dayslog = log(duration_days)

gen countries_deployeddum = 1 if countries_deployed >0
replace countries_deployeddum = 0 if countries_deployeddum  == .

********************************* 表3  多层混合效应负二项回归模型结果（平民伤亡人数）******** 

menbreg deaths_civilians l.deaths_civilians l.countries_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store m1
outreg using "tab5.doc",  replace se var starlevels(10 5 1) sigsymbols(*,**,***)
menbreg deaths_civilians l.deaths_civilians l.units_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store m2
outreg using "tab5.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)
menbreg deaths_civilians l.deaths_civilians l.pko_deployed   l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store m3
outreg using "tab5.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)
menbreg deaths_civilians l.deaths_civilians lsum_notroops  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store m4
outreg using "tab5.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)
************************* 

***************************** 图2  变量的边际效应

menbreg deaths_civilians l.deaths_civilians l.countries_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:

	
gen est = .
gen up = .
gen lw = . 

gen n = 0 if _n== 1
replace n= n[_n-1] +1 if n==.
replace n = . if n>10
	


local j = 1
forvalues i = 0(1)10{
	quietly lincom  l.countries_deployed*`i'  
	replace est = exp(r(estimate)) if _n==`j'
	replace lw = exp(r(estimate)-1.96*.1436874) if _n==`j'
	replace up = exp(r(estimate)+1.96*.1436874) if _n==`j'
	local j = `j' + 1

}

*** 图2-a
graph twoway line lw n, clpattern(dash) clwidth(medium) clcolor(midblue) ///
	|| line up n, clpattern(dash) clwidth(medium) clcolor(midblue) ///			
	|| line est n, color(black) clwidth(medthick)  legend(off) clpattern(solid) ///
	plotregion(fcolor(white)) graphregion(fcolor(white)) ///
	xscale(titlegap(1.5) range(0 10) noext ) ///
	ylabel(0 100 200 300 400 500, nogrid  angle(horizontal) axis(1) labsize(2.5) glcolor(white)) ///
	xlabel(, labsize(2.5) glcolor(white)) ///
    ytitle("平均预测死亡平民人数", size(3) axis(1)) ///
	xtitle("网格部署国家数量", size(3))	///
	title("", position(12) size(3)) ///
	scheme(s1mono) graphregion(fcolor(white) )  plotregion(lstyle(none) margin(l = 0 b = 1))
graph export prediction.png, width(800) height(600)replace	



menbreg deaths_civilians l.deaths_civilians l.units_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:

drop est up lw n
gen est = .
gen up = .
gen lw = . 

gen n = 0 if _n== 1
replace n= n[_n-1] +1 if n==.
replace n = . if n>23
	


local j = 1
forvalues i = 0(1)23{
	quietly lincom  l.units_deployed*`i'  
	replace est = exp(r(estimate)) if _n==`j'
	replace lw = exp(r(estimate)-1.96* .0962558) if _n==`j'
	replace up = exp(r(estimate)+1.96* .0962558) if _n==`j'
	local j = `j' + 1

}

*** 图2-b
graph twoway line lw n, clpattern(dash) clwidth(medium) clcolor(midblue) ///
	|| line up n, clpattern(dash) clwidth(medium) clcolor(midblue) ///			
	|| line est n, color(black) clwidth(medthick)  legend(off) clpattern(solid) ///
	plotregion(fcolor(white)) graphregion(fcolor(white)) ///
	ylabel(0 1000 2000 3000 4000 5000, nogrid  angle(horizontal) axis(1) labsize(2.5) glcolor(white)) ///
	xlabel(, labsize(2.5) glcolor(white)) ///
	xscale(titlegap(1.5) range(0 23) noext ) ///
    ytitle("平均预测死亡平民人数", size(3) axis(1)) ///
	xtitle("网格部署部队数量", size(3))	///
	title("", position(12) size(3)) ///
	scheme(s1mono) graphregion(fcolor(white) )  plotregion(lstyle(none) margin(l = 0 b = 1))
graph export predictionunit.png, width(800) height(600)replace	
************************* ************************* 


*************************  表4  多层混合效应负二项回归模型结果（暴力事件次数）
menbreg violences l.violences l.countries_deployed   l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store m13
outreg using "tab8.doc",  replace  se var starlevels(10 5 1) sigsymbols(*,**,***)
menbreg violences l.violences l.units_deployed   l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store m14
outreg using "tab8.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)
menbreg violences l.violences l.pko_deployed   l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store m15
outreg using "tab8.doc",  merge  se var starlevels(10 5 1) sigsymbols(*,**,***)
menbreg violences l.violences lsum_notroops  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store m16
outreg using "tab8.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)
************************* ************************* 


************************* 表5  双变量Probit模型结果 ************************* 
**** bivariate model
gen deaths_civiliansdum = 1 if deaths_civilians >3
replace deaths_civiliansdum = 0 if deaths_civiliansdum ==.
gen countries_deployeddum = 1 if countries_deployed > 4
replace  countries_deployeddum = 0 if countries_deployeddum ==.

biprobit (countries_deployeddum= instrument pko_africa_total capdist_log  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog) (deaths_civiliansdum=countries_deployeddum  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog ), nolog cl(prio_grid)
outreg using "tabbiviate.doc", se var starlevels(10 5 1) sigsymbols(*,**,***) note (Robust standard errors in parentheses clustered on cell.) replace
gen violences_dum = 1 if violences >0
replace violences_dum = 0 if violences_dum ==.
biprobit (countries_deployeddum= instrument pko_africa_total capdist_log   l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog) (violences_dum =countries_deployeddum  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog ), nolog cl(prio_grid)
outreg using "tabbiviate.doc", se var starlevels(10 5 1) sigsymbols(*,**,***) note (Robust standard errors in parentheses clustered on cell.) merge



*********** 图3  稳健性检验回归系数

label var countries_deployed 网格部署国家数量
label var units_deployed 网格部署部队数量
label var pko_deployed 网格维和平均人数
label var lsum_notroops 网格维和部队总人数（对数）


menbreg deaths_civilians_1 l.deaths_civilians_1 l.countries_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store t1
outreg using "tab_type.doc",  replace se var starlevels(10 5 1) sigsymbols(*,**,***)
menbreg deaths_civilians_2 l.deaths_civilians_2 l.countries_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store t2
outreg using "tab_type.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)
menbreg deaths_civilians_3 l.deaths_civilians_3 l.countries_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store t3
outreg using "tab_type.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)
* binary DV
gen deaths_civiliansdum = 1 if deaths_civilians >0
replace deaths_civiliansdum = 0 if deaths_civiliansdum ==.
melogit deaths_civiliansdum l.deaths_civiliansdum l.countries_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store l1
outreg using "logit.doc",  replace se var starlevels(10 5 1) sigsymbols(*,**,***)

melogit deaths_civiliansdum l.deaths_civiliansdum l.units_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store l2
outreg using "logit.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)

melogit deaths_civiliansdum l.deaths_civiliansdum l.pko_deployed  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store l3
outreg using "logit.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)

melogit deaths_civiliansdum l.deaths_civiliansdum l.lsum_notroops  l.death_totallog l.duration_dayslog  l.pop_log_ l.mnt_ l.ttime_log_  l.nlights_mean l.polity2 l.gdppclog l.poplog || gwno: || year:
estimates store l4
outreg using "logit.doc",  merge se var starlevels(10 5 1) sigsymbols(*,**,***)

*********** 图3-a

coefplot ///
 (t1,  keep(L.countries_deployed)  lpatt(solid) lcol(black)  msym(o) msize(2) mcol(black ) ciopts(lpatt(solid) lcol(black))  mlabposition(1) mlabsize(4.2) mlabgap(*2)) ///
  (t2,  keep(L.countries_deployed )  lpatt(solid) lcol(black) msym(triangle) msize(2) mcol(red ) ciopts(lpatt(solid) lcol(red))  mlabposition(1) mlabsize(4.2) mlabgap(*2)) ///
     (t3, keep(L.countries_deployed ) lpatt(solid) lcol(black)  msym(diamond) msize(2) mcol(gray ) ciopts(lpatt(solid) lcol(gray))  mlabposition(1) mlabsize(4.2) mlabgap(*2)), ///
   scheme(s1mono)  plotregion(lcolor(none))  ///
 horizontal xline(0,lwidth(thin) lpattern(dash)) ///
	grid(none) ///
	levels(95) ///
	legend(size(small) order(2 "模型1：政府-非政府组织之间暴力平民死亡人数" 4 "模型2：非政府-非政府组织之间暴力平民死亡人数" 6 "模型3：单边暴力平民死亡人数") rows(3) region(lcolor(none))) ///
	title("混合效应负二项模型", size(medium)) ///
	ylabel(, notick labsize(medium)) 	
graph export type.png, width(800) height(600)replace	

*********** 图3-b
coefplot ///
 (l1,  keep(L.countries_deployed)  lpatt(solid) lcol(black)  msym(o) msize(2) mcol(black ) ciopts(lpatt(solid) lcol(black))  mlabposition(1) mlabsize(4.2) mlabgap(*2)) ///
  (l2,  keep(L.units_deployed )  lpatt(solid) lcol(black) msym(triangle) msize(2) mcol(red ) ciopts(lpatt(solid) lcol(red))  mlabposition(1) mlabsize(4.2) mlabgap(*2)) ///
   (l3, keep(L.pko_deployed ) lpatt(solid) lcol(black) msym(square) msize(2) mcol(blue ) ciopts(lpatt(solid) lcol(black))  mlabposition(1) mlabsize(4.2) mlabgap(*2)) ///
     (l4, keep(L.lsum_notroops ) lpatt(solid) lcol(black)  msym(diamond) msize(2) mcol(gray ) ciopts(lpatt(solid) lcol(gray))  mlabposition(1) mlabsize(4.2) mlabgap(*2)), ///
   scheme(s1mono)  plotregion(lcolor(none))  ///
 horizontal xline(0,lwidth(thin) lpattern(dash)) ///
	grid(none) ///
	levels(95) ///
	legend(size(small) order(2 "模型1" 4 "模型2" 6 "模型3" 8 "模型4" ) rows(1) region(lcolor(none))) ///
	title("", size(medium)) ///
	ylabel(, notick labsize(medium)) 
graph export logit.png, width(800) height(600)replace	

***************************************** end of code*****************************************


